home *** CD-ROM | disk | FTP | other *** search
/ Aminet 31 / Aminet 31 (1999)(Schatztruhe)[!][Jun 1999].iso / Aminet / dev / c / random.lha / random / random-c / random.c < prev    next >
C/C++ Source or Header  |  1999-02-28  |  5KB  |  325 lines
  1. /*
  2.  
  3. C This random number generator originally appeared in "Toward a Universal 
  4.  
  5. C Random Number Generator" by George Marsaglia and Arif Zaman. 
  6.  
  7. C Florida State University Report: FSU-SCRI-87-50 (1987)
  8.  
  9. C 
  10.  
  11. C It was later modified by F. James and published in "A Review of Pseudo-
  12.  
  13. C random Number Generators" 
  14.  
  15. C 
  16.  
  17. C THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
  18.  
  19. C       (However, a newly discovered technique can yield 
  20.  
  21. C         a period of 10^600. But that is still in the development stage.)
  22.  
  23. C
  24.  
  25. C It passes ALL of the tests for random number generators and has a period 
  26.  
  27. C   of 2^144, is completely portable (gives bit identical results on all 
  28.  
  29. C   machines with at least 24-bit mantissas in the floating point 
  30.  
  31. C   representation). 
  32.  
  33. C 
  34.  
  35. C The algorithm is a combination of a Fibonacci sequence (with lags of 97
  36.  
  37. C   and 33, and operation "subtraction plus one, modulo one") and an 
  38.  
  39. C   "arithmetic sequence" (using subtraction).
  40.  
  41. C======================================================================== 
  42.  
  43. This C language version was written by Jim Butler, and was based on a
  44.  
  45. FORTRAN program posted by David LaSalle of Florida State University.
  46.  
  47. */
  48.  
  49.  
  50.  
  51. #include <stdio.h>
  52.  
  53. #include <stdlib.h>
  54.  
  55.  
  56.  
  57. #define TRUE -1
  58.  
  59. #define FALSE 0
  60.  
  61. #define boolean int
  62.  
  63.  
  64.  
  65. static void rmarin(int ij, int kl);
  66.  
  67. void ranmar(float rvec[], int len);
  68.  
  69.  
  70.  
  71. void main()
  72.  
  73. {
  74.  
  75.     float temp[101];
  76.  
  77.     int ij, kl, len, i;
  78.  
  79.  
  80.  
  81.     /* random seeds for the test case: */
  82.  
  83.     ij = 1802;
  84.  
  85.     kl = 9373;
  86.  
  87.  
  88.  
  89.     /* do the initialization */
  90.  
  91.     rmarin(ij,kl);
  92.  
  93.     
  94.  
  95.     /* generate 20,000 random numbers */
  96.  
  97.     len = 100;
  98.  
  99.     for (i=1; i<=200; i++)
  100.  
  101.         ranmar(temp, len);
  102.  
  103.         
  104.  
  105. /*
  106.  
  107. C If the random number generator is working properly, the next six random
  108.  
  109. C numbers should be:
  110.  
  111. C           6533892.0  14220222.0  7275067.0
  112.  
  113. C           6172232.0  8354498.0   10633180.0
  114.  
  115. */
  116.  
  117.      len = 6;
  118.  
  119.      ranmar(temp,len);
  120.  
  121.      
  122.  
  123.      for (i=1; i<=6; i++)
  124.  
  125.          printf("%12.1f ", 4096.0*4096.0*temp[i]);
  126.  
  127. }
  128.  
  129.  
  130.  
  131. static float u[98], c, cd, cm;
  132.  
  133. static int i97, j97;
  134.  
  135. static boolean test = FALSE;
  136.  
  137.  
  138.  
  139. static void rmarin(ij,kl)
  140.  
  141. int ij, kl;
  142.  
  143. {
  144.  
  145. /*
  146.  
  147. C This is the initialization routine for the random number generator RANMAR()
  148.  
  149. C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
  150.  
  151. C                                                      0 <= KL <= 30081
  152.  
  153. C The random number sequences created by these two seeds are of sufficient 
  154.  
  155. C length to complete an entire calculation with. For example, if sveral 
  156.  
  157. C different groups are working on different parts of the same calculation,
  158.  
  159. C each group could be assigned its own IJ seed. This would leave each group
  160.  
  161. C with 30000 choices for the second seed. That is to say, this random 
  162.  
  163. C number generator can create 900 million different subsequences -- with 
  164.  
  165. C each subsequence having a length of approximately 10^30.
  166.  
  167. C 
  168.  
  169. C Use IJ = 1802 & KL = 9373 to test the random number generator. The
  170.  
  171. C subroutine RANMAR should be used to generate 20000 random numbers.
  172.  
  173. C Then display the next six random numbers generated multiplied by 4096*4096
  174.  
  175. C If the random number generator is working properly, the random numbers
  176.  
  177. C should be:
  178.  
  179. C           6533892.0  14220222.0  7275067.0
  180.  
  181. C           6172232.0  8354498.0   10633180.0
  182.  
  183. */
  184.  
  185.     int i, j, k, l, ii, jj, m;
  186.  
  187.     float s, t;
  188.  
  189.     
  190.  
  191.     if (ij<0 || ij>31328 || kl<0 || kl>30081) {
  192.  
  193.         puts("The first random number seed must have a value between 0 and 31328.");
  194.  
  195.         puts("The second seed must have a value between 0 and 30081.");
  196.  
  197.         exit(1);
  198.  
  199.     }
  200.  
  201.     
  202.  
  203.     i = (ij/177)%177 + 2;
  204.  
  205.     j = ij%177 + 2;
  206.  
  207.     k = (kl/169)%178 + 1;
  208.  
  209.     l = kl%169;
  210.  
  211.     
  212.  
  213.     for (ii=1; ii<=97; ii++) {
  214.  
  215.         s = 0.0;
  216.  
  217.         t = 0.5;
  218.  
  219.         for (jj=1; jj<=24; jj++) {
  220.  
  221.             m = (((i*j)%179)*k) % 179;
  222.  
  223.             i = j;
  224.  
  225.             j = k;
  226.  
  227.             k = m;
  228.  
  229.             l = (53*l + 1) % 169;
  230.  
  231.             if ((l*m)%64 >= 32) s += t;
  232.  
  233.             t *= 0.5;
  234.  
  235.         }
  236.  
  237.         u[ii] = s;
  238.  
  239.     }
  240.  
  241.     
  242.  
  243.     c = 362436.0 / 16777216.0;
  244.  
  245.     cd = 7654321.0 / 16777216.0;
  246.  
  247.     cm = 16777213.0 / 16777216.0;
  248.  
  249.     
  250.  
  251.     i97 = 97;
  252.  
  253.     j97 = 33;
  254.  
  255.     
  256.  
  257.     test = TRUE;
  258.  
  259. }
  260.  
  261.  
  262.  
  263. void ranmar(rvec, len)
  264.  
  265. float rvec[];    /* len random numbers are placed in rvec[1..len] */
  266.  
  267. int len;
  268.  
  269. /*
  270.  
  271. C This is the random number generator proposed by George Marsaglia in 
  272.  
  273. C Florida State University Report: FSU-SCRI-87-50
  274.  
  275. C It was slightly modified by F. James to produce an array of pseudorandom
  276.  
  277. C numbers.
  278.  
  279. */
  280.  
  281. {
  282.  
  283.     int ivec;
  284.  
  285.     float uni;
  286.  
  287.     
  288.  
  289.     if (test==FALSE) {
  290.  
  291.         puts("Call the init routine rmarin() before calling ranmar().");
  292.  
  293.         exit(2);
  294.  
  295.     }
  296.  
  297.     for (ivec=1; ivec<=len; ivec++) {
  298.  
  299.         uni = u[i97] - u[j97];
  300.  
  301.         if (uni < 0.0) uni += 1.0;
  302.  
  303.         u[i97] = uni;
  304.  
  305.         i97--;
  306.  
  307.         if (i97==0) i97 = 97;
  308.  
  309.         j97--;
  310.  
  311.         if (j97==0) j97 = 97;
  312.  
  313.         c -= cd;
  314.  
  315.         if (c<0.0) c += cm;
  316.  
  317.         uni -= c;
  318.  
  319.         if (uni<0.0) uni += 1.0;
  320.  
  321.         rvec[ivec] = uni;
  322.  
  323.     }
  324.  
  325. }